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ABSTRACT 


It is relatively easy to generate, by digital computer, 
large numbers of seemingly independent random numbers with a 
uniform distribution over a fixed range, say ■ | < X n < j-. 

Methods of generating gaussian, or normal, random numbers generally 
are based on either non-linear transformations on random numbers 
from a uniform population, or the summing of enough independent 
numbers from a uniform population for the central limit theorem 
to be applicable. In the first case a time-consuming evaluation 
of a complicated function is involved. The second method is also 
slow because a large number of uniform random variables must be 
generated and summed for each normal random variable obtained. 

This note discloses a method based on the central limit theorem, 
except that the summing of N uniform random variables gives N 
normal random variables. The approach is to form an N dimensional 
vector whose components are uniform random variables, multiply the 
vector by a Hadamard matrix, and use the resulting components as 
normal random variables. It can be shown that the resulting N 
components have a uniform density inside a N-dimensional hypercube 
aligned with diagonals along the coordinate axes. However, the 
one dimensional marginal densities, the two dimensional marginal 
densities, indeed all the marginal densities tend toward the normal 
density as N gets large. Furthermore the components are uncorre¬ 
lated and have equal variance, independent of N. However, some 
of the fourth moments, which should be zero for independent 
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normal random variables, are not zero for our derived set 
(although these moments do approach zero as N becomes large). 
These moments can be made zero, however, by randomly changing, 
or not changing, the sign of each component. 

The method proposed is very fast because the principal 
step, Hadamard matrix multiplication, requires only N log2 N 
additions to produce N components. 


Accepted for the Air Force 
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A NEW METHOD OF GENERATING GAUSSIAN RANDOM 
VARIABLES BY COMPUTER 

Many scientific computations require the generation by a 
computer of a large number of seemingly random numbers^ with a 
multidimensional normal (gaussian) probability density function. 
However, although there are efficient algorithms available for 
generating uniformly distributed random numbers, the algorithms 
which generate normal random numbers are generally much slower. 

One typical algorithm begins by generating a moderate number of 
uniform random variables (say N) and forming the sum, which, for 
large N, is approximately normal, according to the central limit 
theorem. This method is slow because many uniformly distributed 
random variables must be developed for each resulting normal 
random variable, and because many additions must be performed 
for each normal random variable, as well. Another typical method, 
more attractive because it is less dependent on approximation, 
begins with a single uniform random variable and performs a non¬ 
linear transformation of the random variable to form a new random 
variable with a normal density. Unfortunately, the non-linear 
function which must be evaluated is complicated and cannot be 
computed quickly on most digital computers. 

This note describes a variation of the central limit theorem 
approach, in which N uniform random variables (r.v.) are trans¬ 
formed into N approximately normal r.v. by a Hadamard transformation. 

rr The term "pseudorandom" is often used to describe the determin¬ 
istically generated, but seemingly random numbers. 
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This operation takes ^ N log 2 N additions and ^ N log 2 N subtractions, 
which is to say, log 2 N additions or subtractions per normal r.v. 
obtained, and, of course, it is only necessary to generate one 
uniform r.v. for each normal r.v. obtained. The plan of the 
paper is to first describe the ideal uniform and normal multivariate 
densities, second to describe the Hadamard transformation used, 
third to derive the properties of the transformed variables, and 
fourth, to consider the properties of the transformed variables 
after they have been subjected to random sign changes. 


I. Properties of the Ideal Uniform and Normal Densities 


We define the ideal uniform r.v., X n , as one whose probability 


density is 


w - 


1 

0 


xj < 


n 


a) 


IV > 7 

For several such r.v.s Xq,X^,...,X N _^ to be independent implies 
that the joint density is 

1 if all |X n l < \ 


(X n ,Xi ,x 9 ,... >x^_-^) 


xx...x v 0* 1* 2 


0 if any |X n l > j 


( 2 ) 


On an N dimensional space, the joint density is unity on the 
inside and zero on the outside of a hypercube centered on the 
origin with coordinate axes passing from the center of each face 
to the center of the opposite face. 
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The expected value, E(X n ) = X n , 

all the other moments can be found. 

moments of X , 
n 7 


is zero and, by integration, 
For the case of even order 


= 1 

n 2 2p (2p+l) 


and for odd order moments 


(3) 


X 2p+1 = 0 


(4) 


Since the mean is zero, the concept of central moments is super¬ 
fluous. For joint moments, we appeal to independence, which says 
that the moment of a product is the product of the moments. 
Therefore 


(X a X b X c ) = 0£>( x S> (X c> (5) 

and, of course, these joint moments are zero except when m, n, Z 
are all even. 

It is worth taking note of the second moment, X n = 1/12. 

The normal probability density of an r.v. Y , with zero mean, 


is given by 








2 

where a is the variance. We shall be comparing several densities 
in this note and we shall force the variances to all be equal to 
1/12. Therefore a 2 = 1/12. 

The multivariate density of several such independent variables, 
Y 0’ Y 1>••• ,Y N-1> ls 


P yyy--y (v o’ Y i’"' ,Y N-i ) = 


f7v\ 


1 -( Y o +Y l + "- +Y N-l) /2 ° 2 

n fT e 


tt a 


(7) 


The odd order one dimensional moments of Y are 

n 


Y 2p+1 = 0 
n 


( 8 ) 


while the even order moments may be shown to be 


p = 1*3 *5•••(2p-l) = l-3»5’••(2p-l) 


(1/cr) 


2 ? 


12 p 


(9) 


Again, by independence, the joint moments are equal to the product 
of the individual moments, so that the moment 


Y X^ = o£)0£)( Y i) 


a y v b y v c 


( 10 ) 


is zero unless m, n, SL are all even. 
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The normal density is never zero, and is, in fact, radially 
symmetrical, i.e., P y(Y q,Y^, . . . ,Y^_-^) '*' S on -*-y a f unct ion of 

the distance of Yq,Y^, . . . from the origin in N-space. This 

is also a true statement for any marginal density. The same 
statement is not true for the multivariate uniform density, which 
takes on different probabilities in the directions of the vertices 
of the hypercube than along the direction of the faces of the cube. 

II. The Hadamard Matrix 

By definition, a Hadamard matrix is an orthogonal matrix 
whose elements are all either +1 or -1. Hadamard matrices do not 
exist for all dimensionalities; for example there is no 3 x 3 
Hadamard matrix. It is relatively easy to construct an N X N 
Hadamard matrix, however, if N is a power of two. We shall content 
ourselves with a particular Hadamard matrix for each power of 
two, for our purpose is not to investigate Hadamard matrices but 
to apply them. Letting N = 2 , and numbering our dimensions in 
the binary system 

i = i~ + 2i, + 4i 0 + ... + 2 c “ 1 i , 0 ^ i s N-l 

(J 1 z c-1 

( 11 ) 

j - j 0 + + 4j 2 + ... + 2 C ‘ 1 j c . 1 0 i j i N-l 

the elements of our particular Hadamard matrix, H, are given by 
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h 1 , = (-l) 1()Jo (-l) llJl (-l) l2J2 

X 5 J 


(-D 


1 c-l-*c-l 


( 12 ) 


We shall show in an appendix that multiplication of an N component 
vector by H can be accomplished with N log 2 N additions and 
subtractions. 

It is clear that h. . = h. . so that H is a symmetric matrix. 

1 > J J > 1 

ttl ttl 

Also, the zero c row and the zero c column are composed of +l's 

j. ^ j . 

only. The product of the elements in the i c row with the k c row 
give the elements of another row of the Hadamard matrix. 



k,j 


= h 




where i is the row whose number is 


(13) 


l = (i 0 © k Q ) + 2(i L © k x ) + ... 

+ 2 c - 1 (i c . 1 9 k c . L ) (14) 

where © means "exclusive or". We shall say as a shorthand, 

l = i © k 

in place of (14). By symmetry, the product of columns is also a 
column. 
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A consequence of (13) is the orthogonality of the Hadamard 


matrix. Strictly speaking the matrix 


1 

7n 


H 


is not only its own transpose but its own inverse, and it is the 
matrix which we shall use in the transformation of uniform r.v. 
to normal r.v. 

We give as an example the 4x4 Hadamard matrix 


H 


1 

1 

1 

1 


1 1 

-1 1 

1 -1 

-1 -1 


1 

-1 


-1 

1 


III. Hadamard Generation of Almost Normal r.v. 


We begin with a set of uniformly distributed r.v. assumed 
independent. We form a new set of r.v. from the uniform set by 
the operation 


W — -ttt 

m /N 


N-l 
Z h 
n=0 


X 


n,m n 


m = 0,1, 


.N-l 


(15) 


and we shall be interested in the properties of the set of r.v. 

W . For any single random variable, it is easy to see that P w (W m ) 
is the density of a sum of independent r.v. with zero mean and 
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equal variance. It is known from the central limit theorem that 
such a sum approaches a gaussian r.v. in the limit of large N. 


The approximation is good in the region of W m near zero, bad in 
the tails. It is, of course, clear that 



(16) 


so that P w (W m ) is zero outside this range. However, for large N, 
2 /N is many standard deviations from the mean, so that the 
difference can be tolerated. Indeed, the approximation is the 
same as for the common method of generating normal r.v. except 
that we hope to use a larger value of N and so obtain a good 
approximation. 


We have generated, however, N different r.v. from the same 


N uniform r.v. and we should not expect the W m 's to be independent. 
This leads us to investigate the multivariate density 

P ww.. .w^ W 0 ,W l’ * * * » W N-1^ * But 

( w 0 ,W l» * ’ ’ = P xx...x (X 0’ X l’-"’ X N-l)/ 


ww, 


,w 


where 


D 


BW o 

dX o 



£ 

N-l 


aw. 

ax 


N-l ,..., 
0 


aw 

ax 


N-l 

N-l 


p 

ww...w 


(w 0 ,w 1! 


,W N-1^ 


P 

XX... X 




• • 


• }X N-1 


) 


(17) 
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since the magnitude of the determinant of H is unity. Of 
course, in (17), P (X n ,X,,...,X M ,) is to be evaluated at 

XX • • • X U X 1 N*“J- 


' x o ' 
X 1 

1 ,T 

_ w o ' 

w i 

• 

" Tfi H 

# 

(- 

X 

> 

.j 

i _ 


1 

• * 

> 

mA 

1 _ 


but this is only a change of coordinate system. In fact, there¬ 
fore, the N-dimensional joint density P is not normal but 

uniform. The density is one inside and zero outside a hypercube, 
This hypercube, however, is rotated. Whereas for the variables 
X , the hypercube was situated with the coordinate axes passing 
through its faces, for the W n 's we shall see that the coordinate 
axes pass through vertices of the hypercube. In fact, the 
coordinates of some of the vertices of the hypercube surrounding 

P xx...x^ X 0 ,X l»••• ,X N-1^ are 


Xj = y 


h. . 


for any i, and such a vertex becomes rotated into the position 
with coordinates 


W. 

J 


| yrr 


0 


j = i 
j * i 
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ttl 

which lies on the coordinate axis of the i c coordinate in W-space. 
The vertex at -Xj rotates to the opposite position, thus proving 
that the coordinate axes of W-space pass through diagonals of the 
hypercub^. 

A hypercube in N-dimensions has 2^ vertices. There are, 
however, only N coordinate axes, so only 2N vertices can lie on 
the coordinate axes. The remaining 2^ - 2N vertices of the 
hypercube do not map onto the coordinate axes. At least some of 
the vertices are found in the "interior" of hyperquadrants, but 
not every hyperquadrant of the N-space can be lucky enough to have 
a vertex within it, for there are 2^ hyperquadrants and at most 
2^ - 2N vertices to go around. Therefore, we learn that 
P^ w (Wq,W^,...»^n_i) is not only radially asymmetrical but is 
not symmetrical in each hyperquadrant of W-space. This shows up 

when we compute the moments of P . However, the most important 

ww..,w * 

moment of P , namely W W , is zero, showing that the variables 

ww..,w* 3 m n* * ° 

W are at least uncorrelated. We shall now compute all the first 
m 

and second moments of P w , (W n ,W-, , . . . ,W M -, ) . 

ww...w U* 1 1 N-l 

The first moments are all equal and are all zero. 


W = -7TT 

m /N 


E h 

r 

n 


X n - 7N 


E h 

I 

n 


X n = 0 


( 18 ) 
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because X =0. The second moments are 
n 


N-l N-l 


W W, = ± E Eh h 
ah N _ n r. m 
m=0 n=0 


,a n,b Vn 


(19) 


but X m X n is zero unless m = n, in which case it is a . Thus (19) 
becomes 


_ 1 N-l 9 

W W, = rr E h / \ a 

a b N m, (a©b) 

m=0 


( 20 ) 


th 


By definition of a © b, it is zero if a = b, and the zero column 
of a Hadamard matrix consists of N ones. Therefore 


r 7i 2 

W = a 

a 


( 21 ) 


as expected. But if a / b, a 0 b is not zero and sirce all the 
other columns of H contain as many -l's as l's, the sum 


E h. 

m 


l m, (a©b) 


is zero, giving 


W W, 
a b 


= 0 


as stated above. 


( 22 ) 
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It is not difficult to argue that all the odd order moments 
of P rrrT __ are zero. We focus now on fourth moments. By definition 

WW...W J 


W W K W W , = —nr 
abed n 2 


h. h. , h 

j ,b 


N-1 

N-1 

N-1 

N-1 

X 

X 

X 

X 

i=0 

j=0 

O 

II 

X 

1=0 


, h . , X.X.X, X. 

k,c £,d 1 j k l 


(23) 


To evaluate the fourth moments we note that X.X.X, X. = 0 for 

i j k a 

most combinations of i,j,k,£ with the exceptions listed below: 

1 


1. 

2 . 

3. 

4. 


T7I 


a = 


144 


i=j , k=j i, ij^k for which X^X^ = 

i=k, j = £, i^j for which X^xf = cr = 

i=£, j =k, i/j for which X?X? = 


i=j=k=£ 


~~K 

for which X. 


1 

80 


These are all disjoint cases. Therefore 


W 


1 7 1 


2/n-I 


/N-1 


a W b W c W d = ^ [ ± l Q h i,a0b R f 0 h k,c0d] 

Vk^i 

i i 2 /n-I /n-1 \ 

+ 7 (T?> li=o hi . a& (j=o hk . b0d 

X Vi ' 


1 7 1 


2/n-i 


'n-1 


+ 7 (t?) U-0 hi > a ® d ij=0 h j > b 6cJ 

Vi 


N-1 


+ ^7 ± l 0 h i,a0bfc0d 


(24) 
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Equation (24) looks very complicated but it may be straightforwardly 

~T 

evaluated for any special case. Thus, for W we have a=b=c=d and 

ci 

we find 


W? - ~T (A) (N 2 -N) + ~T (tj) (N 2 -N) + 

a N / iz N z iz 

(vf) (n 2 -n) + ($7j) N 

N N OU 

^ = WT + ^ ( M “ wr) = 3 ° 4 ' IM (25) 

Similarly, if we let a=c, b=d, a^b we get 

" ^7 <T7> (n 2 - n ) + ^7 ( ST> )N 

- Ttt + Ton - ° 4 + Bn < 26 > 

These results approach the gaussian independent results 3o\ a^ 
for large N. But consider the curious case of a,b,c,d all unequal 
but with a0b0c0d = 0. For this 


W a W b W c W d " SM 


(27) 
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An example of such a fourth moment is WqW^W 2^3 • This 
moment also approaches the expected result for independent normal 
r.v. as N - 00 . However, it is somewhat disturbing that it is 
non-zero. We can, after all, convince ourselves that if a 
certain moment is within, say, 27 o of the expected result for 
normal independent r.v. that this is somehow good enough. But 
what can we say if the expected moment should be zero? It is 
not clear how important this is, and indeed it probably depends 
on the application for which the r.v. are desired. 

Before considering the effect of a random sign change, it 
is well to summarize what we have shown so far. We have proposed 
a method by which uniformly distributed zero-mean random noise 
can be converted to approximately normal zero mean random noise. 
The one dimensional marginal densities are sums of N independent 
uniform r.v. and thus tend nicely toward normal for large N. 

The r.v. are also uncorrelated (in the sense that the expected 
value of a product is zero). Furthermore, based on the moments 
we have computed, the difference between a given joint moment 
and the same joint moment for independent gaussian noise seems 
to go to zero as 1/N. N can reasonably be made quite large since 
the amount of computation is proportional to log 2 N per normal 
r.v. generated. 
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IV. HADAMARD GENERATION OF MORE NEARLY INDEPENDENT RANDOM 

VARIABLES 


Let us now suppose that we have generated N random variables 

W . We have seen that some of the moments W W U W W, , which 
m abed’ 

would be zero for an independent gaussian set, are non-zero. We 
can form a set of N new random variables from the W's by changing 
the sign of each W, or not changing it, at random. That is, 
consider an easily obtained set of r.v., r m , independent of each 
other and of W m , with equal likelihood of being +1 or -1. For 
these r.v. the moments are all either zero or one. Even order 
moments are equal to 1, 



odd order moments are equal to zero, 


2 p + 1 

r r 


= 0 


5 


and, by independence, the moment of the product of different r m 's 
is the product of the moments. Now consider the set of random 
variables 


V = r W 
m mm 


(28) 
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We shall compute some of the moments of the joint density 


P vv v (Vq, ^1 ,,, '*^N-1^ to s ^ ow that this density is nearly 
gaussian, for large N. First consider the moments which 


violated the gaussian assumption for the W's. The is 

equal to 


Wc'dWcV < r a r b r c r d> WcV 


and unless a, b, c, d are equal at least in pairs, r a r b r c r d will 
be zero. A similar argument will show that all joint moments of 

P which would be zero for independent, zero mean gaussian 

noise are indeed zero. However, all other moments will be 

—7 2 

unaffected by the conversion of W to V . Thus V = cr , 

3 m m a ’ 

V V, = 0, and 
a b ’ 


„ 2„ 2 , T 2.. 2 4 , a 

V a V b " W a W b * a + T5N 



( 29 ) 


We next consider the 



order moments 















Discarding the trivially zero moments, we have for the other 


V Vi V V j V V.p 
a b c d e f 


1 

7 


SELSLL 
i j k ( m n 


h. h. , 

i, a J>b 


V 


h . A 

c £,d 


h h . X.X.X, X.X X 
m,e n,f i j k £ m n 


We must again consider moments X.X.X,X.X X , most of which are 

1 j K £ m n ’ 

zero. The exceptions are listed below with the equal subscripts 
grouped and the ungrouped subscripts unequal. Thus 


Disjoint Cases 


Moment 


i j k £ m n 


1 

448 


i 

j 

k 

J l 

m 

n 


i 

j 

k 

m 

£ 

n 


i 

j 

k 

n 

£ 

m 


i 

j 

£ 

m 

k 

n 


i 

j 

£ 

n 

k 

m 

i 








i 

k 

£ 

n 

j 

n 


i 

k 

£ 

n 

j 

m 


j 

k 

£ 

m 

i 

n 


j 

k 

£ 

n 

i 

m 


k 

l 

m 

n 

i 

j 


j 

£ 

m 

n 

* 

i 

k 


j 

k 

m 

n 

i 

j l 


i 

k 

m 

n 

j 

l 


i 

£ 

m 

n 

• 

J 

k 


i 


m 

n 

i 

j i 
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Disjoint Cases 


Moment 


i j 
i j 
i j 

i k 
i k 
i k 
i l 
i l 
i l 
i m 
i m 
i m 
i n 
i n 
i n 


k l 
k m 
k n 
3 l 
3 m 
j n 
j k 
j m 
j n 
j k 

j ^ 
j n 
j k 

j * 
3 m 


m n 
n t 
m l 
m n 
l n 
t m 
m n 
n k 
m k 
t n 
k n 
k l 
t m 
k m 
k l 


1 

1728 


t h 

This table allows us to compute all the 6 C order moments 
for special cases. For example, if we were to compute V ^ , we 

a 

would have as part of the computation a sum 


1 


N-l 

L 

i=0 


h i, 


(a®a@a6a®a 


7^ 

0 a) X i 


1 1 
448 ^ 


and fifteen sums similar to 


18 





I N-l 

^i(a © a © a © a) 


N-l 

S 0 h 

m=U m, 
m^i 


a 


0 a 



1 

c^> 

N J 



■ 95TT 



fifteen sums 

similar to 



, N-l 


N-l 

N-l 

—J s 
N J i=0 

h i, (a © a) 

L h, 
k=0 

(a 0 a) m=0 



k^i 

m^i 




m^k 


2 2 2 

[ i X k X m 


1 

1728 



t h 

Gathering all the terms together, we get, for the 6 C moment 



15 

1728 


N (N-l) (N-2) x 15 /N (N-l) * , 1 ,N, 

-^ + 960 + 448 


V 

a 


15cr + terms 


in r? and 
iN N z 


(30) 


— 5 —7 

Evaluating V V, we set, for example, a = c = d = e, b = f, 

8. D 

a ^ b. Therefore, the sum 


19 













is the same as before, but many of the other thirty sums are 
zero, leaving only the cases 


i 

j 

t 

n 

k 

m 


i 

k 

j 

n 

t 

m 

i 

k 

l 

m 

j 

n 


i 

l 

j 

n 

m 

k 

j 

k 

l 

n 

i 

m 

and 

i 

m 

j 

n 

k 

t 

j 

l 

m 

n 

i 

k 








j 

k 

m 

n 

i 

l 








i 

j 

m 

n 

k 

l 









We quickly count 



3 

1728 


( N(N-l)(N-2) ) 

N 3 


960 


/N (N-l)v 1 ,N v 

( ' 3 ') + TTq C-t) 

N J 440 




+ terms 


• 1 1 

in rr and — 

N ? 


(31) 


— 2 — 2—2 

Similarly for V V u v we obtain 
a b c 


— 2 — 2 — 2 " 

v v u v 

a b c 


— 2 — 2 —1 

V„ V K Z V z 
a b c 


_ ( , N(N-l) (N-2) v 3 /NCN-D v . 1 ( N ) 

" 1728 C > + 95U C } + 


+ terms in rr and -i- 
W N Z 


( 32 ) 















These results also approach the expected results for 

6 6 6 

normal r.v. with an error like 1/N, i.e., 15o , 3cr , and a . 

At this point, we comment that the behavior, for large N, 
of any moment is a matter of combinatorial analysis. For example, 
the moment 


V a P = 1*3*5 • •• (2p-l)+ terms in ^ , etc. (33) 

can be derived by a counting process. For more complicated 
moments, the analysis is so complicated that we have not under¬ 
taken it. 

We previously described the N dimensional distribution of 

the W's in terms of a hypercube which we argued was asymmetically 

N 

distributed among the 2 hyperquadrants in N space. We cannot 
make a simple statement about v (Vq, because 

it is not uniformly distributed, but rather has some density 
which is a continuous function of spacial coordinates. However, 
it is clearly symmetrically distributed in the hyperquadrants, 
by construction. 

Programming Considerations 

The method proposed requires at least the following two or 
three steps: 
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1. Compute N uniform r.v. Xq, X n _^. Suppose that each 

r.v. requires a time T u to generate. Then part 1 will require 
a time NT u * 

2. Compute the "Hadamard transform" of the array of N numbers 

from part 1. This can be done in-place (i.e., without any extra 

memory) with N log 2 N additions or subtractions. Assuming that 

an addition and a subtraction take the same amount of time T , 

a’ 

part 2 will require a time N(log 0 N)T . If it is desired to 

Z 3. 

adjust the variance of the resulting W^, we should add to this 

time an additional time N T , where T is the time for a multiply 

m’ m r J 

(or a scale if that is acceptable). 


3. If the improved approximation to independent gaussian noise 

is desired, generate N additional random variables r^ which can 

be used to randomly change or not change the sign of the r.v. 

Wj generated in part 2. This should take at most a time 

N(T + T ) assuming the r. can be generated as quickly as a 
u a j 

uniform r.v. and sign changing can be accomplished as quickly 
as addition or subtraction. 

The total time, therefore, is like 


T total * N < 2T u + 1 o S2< 2N > V 

but this is the time for generating N r.v. whereas the relevant 
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time would be the time per random variable. This latter time 
is 


Tg - 2 T U + U°g 2 N + 1) T, . 

If many r.v. are needed, as will usually be the case, they 
should be generated N at a time by a subroutine. This sub¬ 
routine could, of course, supply them one at a time to a user 
program, recomputing a block of N when its supply was exhausted. 
Thus, unless fewer than N r.v. were required, the limitation on 
computing the r.v. N at a time should not be important. 

Conclusions 

Random variables with an almost normal distribution can be 
computed by first obtaining uniform random variables, using 
standard computer algorithms, and then performing a linear trans- 
formation on a set of N = 2 of the uniform random variables to 
obtain N random variables whose joint distribution is not 
gaussian but uniform in a hypercubic region. Nevertheless, many 
of the marginal densities are nearly gaussian, in accord with 
the central limit theorem, and the r.v. are uncorrelated. In 
the limit of large N, all the moments which have been determined 
seem to approach the moments expected for gaussian independent 
r.v. However, some of the joint moments which would be zero for 
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independent gaussian r.v. are finite for the derived set. By 
randomly changing or not changing the sign of the linearly trans¬ 
formed r.v., a new set of r.v. is obtained for which all the 
moments, including joint moments, show considerable similarity 
to the moments expected from the gaussian independent case, for 
large N. The procedure has the advantage of being faster than 
other decent methods of generating gaussian r.v. on a digital 
computer, with relatively good accuracy. 
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APPENDIX 


NUMBER OF OPERATIONS IN HADAMARD MATRIX MULTIPLICATION 

If we consider the multiplication of an arbitrary vector by 

the matrix H, with elements h. . given by (12), an abvious 

1 9 J 

proceedure is to partition the matrix into four quadrants, 


H = 


l. • 


c . . 


b. . 


d ij 


a M 

= h. . -s 



b U 

= h. . N 

+ J 


o.i.f- 

C i.j 

= h. , N . 
i + j,J 

> 

0 s j * f - 

d i,j 

= h. , N . , N 

i + 7 ,j + 7 



i that 

for 



a. . 

f (p 0 ^1-^1 

= (-1) 0 0 (-1) 1 1 ... 

(-l) ic -2 J c-2 


(-i) 


: “l^C -1 


i c _l = J c _i = 0 for t ^ ie entire quadrant. Therefore a^ ^ is 


actually an ^ x ^ Hadamard matrix. Similarly 


a • • b • • c• . -d. • 

i,j i.J 
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It follows that if the vector to be multiplied by H in a 
column vector [X^] we may form two shorter vectors with components 


I 

X i 


It 

X i 


X, + X. , N 

1 1 T *2" 


X. - X. , N 
l 1 + 7 


0 ^ i * 


0 * i * 


N 


N 


- 1 


- 1 


N 

and we shall find that the first ^ components of the product are 
given by 

N 


E 

i=0 


a. . 


N 


X, 


and the second components are given by 
N 


T l 

E a. . 

i-0 i *J 


ft 

X. 


N N 


each of which are x Hadamard matrix multiplications. This 
reduction cost us N/2 additions and N/2 subtractions. Of course 


a. . may be similarly partitioned and this may be done repeatedly. 
1 y 3 

N N 

Each reduction costs and additional -jr additions and ^ subtractions 
so that a total of Nc = N log 2 N operations is required in all. 
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